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CONVECTION-DIFFUSION-REACTION EQUATIONS* 

SHAOHONG DUt AND XIAOPING XIE* 

Abstract. We prove the convergence of an adaptive mixed finite element method (AMFEM) for 
(nonsymmetric) convection-diffusion-reaction equations. The convergence result holds from the cases 
where convection or reaction is not present to convection-or reaction-dominated problems. A novel 
technique of analysis is developed without any quasi orthogonality for stress and displacement variables, 
and without marking the oscillation dependent on discrete solutions and data. We show that AMFEM 
is a contraction of the error of the stress and displacement variables plus some quantity. Numerical 
experiments confirm the theoretical results. 
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1. Introduction and main results. Let il be a bounded polygonal or polyhedral domain 
in R d , d = 2 or 3. We consider the following convection-diffusion-reaction equations: 


where £ e L°°(Q;R dxd ) is an inhomogeneous and anisotropic diffusion-dispersion tensor, 
w is a (dominating) velocity field, r is a reaction function, and / is a source term. The choice 
of homogeneous boundary conditions is made for ease of presentation, since similar results 
are valid for other boundary conditions. 

Adaptive methods for the numerical solution of PDEs are now standard tools in science 
and engineering to achieve better accuracy with minimum degrees of freedom. The adaptive 
procedure of (1. 1) consists of loops of the form 


A posteriori error estimation (ESTIMATE) is an essential ingredient of adaptivity, and reaches 
its mature level after two decades of development [1, 3, 4, 5, 1 1, 12, 18, 38]. However, the 
analysis of convergence of the whole algorithm (1.2) is still in its infancy, and is carried out 
mainly for standard adaptive finite element methods (AFEM) [8, 23, 30, 31, 32]. 

Due to the saddle-point characteristic of mixed finite element approximation, there is no 
orthogonality available, as is one of main difficulties in the convergence analysis of AMFEM. 
Thus one has to find some quasi-orthogonality instead of the orthogonality, and the occur- 
rence of oscillation of data is inevitable. Hence, how to deal with the oscillation becomes a 
key issue in the analysis. For the convergence of AMFEM, the present studies mainly focus 
on Poisson equations. In [10], Carstensen and Hoppe proved the error reduction and conver- 
gence for the lowest-order Raviart-Thomas element with marking the data oscillation. Chen, 
Hoist andXu [14] showed the convergence of a quasi-error with marking the data oscillation. 
In [6, 13, 19], the convergence was analyzed for the lowest-order Raviart-Thomas element 
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V • (SVp) + V • (pw) +rp= f in ft, 
p = on <9ft 


(1.1) 


SOLVE -> ESTIMATE MARK REFINE. 


(1.2) 


1 
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where the local refinement was performed by using only either the estimators or the data 
oscillation term. 

For general diffusion problems and more general mixed elements, by using the orthogo- 
nality of the divergence of the flux, Du and Xie [20] showed the convergence of the flux error 
plus some quantity without marking the oscillation. 

The purpose of this paper is to prove the following convergence results for an AMFEM 
for the convection-diffusion-reaction equations (1.1) and verify them computationally. 

Theorem 1.1. (Convergence of AMFEM) Denote by {T k , e\, A 2 k } k > the sequence of 
meshes, the error of the stress and displacement variables, and some quantity (defined by 
(5.12)) produced by the AMFEM algorithm. Let h be the mesh size of the quasi-uniform 
initial mesh %■ Then there exist two positive constants q and a e (0, 1) such that 

4+1 + (1 - h q)A 2 k+1 < a 2 (e 2 k + (1 - h q)A 2 k ) 

when Iiq < *~° 2 K This means that AMFEM, as h a is small enough, converges with a linear 
rate a, namely, 

e 2 + (l-h q)A 2 < a 2k (e 2 + (l-h q)A 2 a ). 

This theorem extends the convergence results in [20] in the following several aspects. 

• We deal with more general convection-diffusion-reaction equations here with vari- 
able coefficients S, w and r, whereas in [20] w and r vanish. 

• The orthogonality for the divergence of the flux is absent due to the convection term 
w • Vp and the zero order term (r + V • w)p. So this contribution considers not only 
the flux (stress variable) error but also the displacement variable error. 

• The quasi-orthogonality for stress and displacement variables also fails due to the 
terms w • Vp and (r + V • w)p. This will lead to an additional constraint on the 
mesh size, ho, of the quasi-uniform initial mesh %■ 

• The oscillation term depends on the discrete solution and data. Therefore, the oscil- 
lation and error can not be reduced separately here. In [20] the oscillation term is 
not included in the a posteriori indicators. 

• Since the error and oscillation are now coupled, in order to prove convergence with- 
out marking the oscillation term, we need to handle them together by following the 
same idea as in [16, 26]. 

• In comparison with previous analysis methods, ours considers the a posteriori indi- 
cators with weighted factors. We also release the constraint that the divergence of 
the convection term is free in contrast to the analysis of standard AFEM (see [27]). 

The rest of this paper is organized as follows. Section 2 gives some preliminaries and 
details on notations. Section 3 derives an estimate for the error between L 2 — projection of 
the displacement and its approximation solution, which is key to the convergence analysis. 
Section 4 shows the estimator reduction. We prove theorem 1.1 (Convergence of AMFEM 
algorithm) in section 5 and present four numerical experiments to illustrate properties of 
AFMEM in section 6. 

2. Assumptions, weak problem, and AMFEM algorithm. For a domain A c R d , 
we denote by L 2 (A) and L 2 (^4) =: (L 2 (A)) d the spaces of square-integrable functions, by 
(•, -)a the L 2 (A) or L 2 (^4) inner product, by || • 1 1 ^ the associated norm, and by \A\ the 
Lebesgue measure of A. Let H k (A) be the usual Sobolev space equipped with norm 1 1 • \ \ k .A 
for k = 1,2; H%(A) := {v e H X (A) : v\ dA = 0}, H(div,A) := {v e L 2 (A) : divv'e 
L 2 (A)}. < ■, ■ >qa denotes d — 1-dimensional inner product on dA for the duality paring 
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between H 1 l 2 {dA) and H 1 ^ 2 (dA). In what follows we shall omit the subscript f2 when 

A = n. 

Let Th be a shape regular triangulation in the sense of [15], and denote the mesh size 
h T := |T| 1/,d with |T| the area of T £ Th- Let Cq be a positive constant depending only on a 
quantity Q, and C{(i = 1, 2, • • • ) positive constants determined only by the shape regularity 
of Th- We denote by £/, the set of element sides in 7/», by e° the set of interior sides of 
elements. For K e Th, denote by sk the set of sides of K. Furthermore, we denote by 
ujk and u a the unions of all elements in Th respectively sharing a side with K and sharing 
a side a e e h . We use the "broken Sobolev space" H 1 (\JTh) ■= {<p € L 2 (Q) : ip\ K e 
H^K^VKeTh}. H 2 ({JT h ) is denned analogously. Denote by [v]\„ := - 
the jump of v € i? 1 (U 7/») over an interior side a := K n L of diameter := diam(cr), 
shared by the two neighboring (closed) elements K, L e Th- Especially, [v]\ a := (u|k-)|<t if 
cr C dK n c?Sl. Note that [•] is a linear operator over the broken Sobolev space -ff 1 (U Th)- 

We note that throughout the paper, the local version of differential operator V is un- 
derstood in the distribution sense, namely, : -ff 1 (U7ft) ~~ ^ {L 2 {&)) d is defined with 
V h v\ K := V(v\ K ) for all KeT h . 

Given a unit normal vector n CT = (m, • • • , rid) T along the side a with d = 2, 3, we 
define the tangential component of a vector veR 11 with respect to n CT by 


where x denotes the usual vector product of two vectors in R 3 . 

Following [39], we suppose that there exists an original triangulation % of ft such that 
data of the problem (1.1) are given in the following way. 
Assumptions of data : 

(Dl) Sk '■= S\k is a constant, symmetric, and uniformly positive definite tensor such that 
cs,kv • v < Skv ■ v < Cs.kv ■ v holds for all v e R d and all K e % with cs,k > 
0, Cs,k > 0; 

(D2) w e RT (To) (see below) and \w\k\ < C Wi k for all K £ To with C^.k > 0; 
(D3) tk '■= r\n is a constant for all K € To; 

(D4) Cw^k ■= 1/2V • w\k + r\ K > and C^-^k ■= |V • w\k + r K \for all K £ To; 
(D5) f e L 2 (Q); 

(D6) ifew.r^K = 0, then C Wirj /( = 0. 

Note that in [21, 22] J\k is assumed to be a polynomial of degree at most k for each 
K e To so as to derive the efficiency of the residual indicators. Here we relax the restriction 
of/(cf. (D5)). 

Introduce the stress variable u := — SVp, the mixed variatinal problem of (1.1) reads as: 
Find (u,p) e i?(div,Q) x L 2 (il) such that 


(V -u,ip)- (S^u ■ w, ip) + ((r + V • w)p, v?) = (/, for all p e L 2 («). (2.2) 

Let Po(K) denote the set of constant functions on each K e Th- We respectively define 
the lowest order Raviart-Thomas finite element ([35]) space and the piecewise constant space 
as following: 



(S^u, v) - (p, V ■ v) = for all v e (div, O), 


(2.1) 


q h G H(div, O) : Vif G Th, 3a e K d , 3& G K, 
such that q/i(x) = a + 6x, for all x e if. 
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and 

Po(%) := {v h G : G %, v h \ K G P (*T)}- 

We note that V • (RT (T h )) C P (T ft ). 

The centered mixed finite element scheme (cf. [17, 39]) of (l.l)reads as: Find (uh,Ph) € 
PT (T fe ) x P (% ) such that 

{S- 1 u h ,v h )-(jp h ,V-v h ) = forallv ft Gi?T (7; i ), (2.3) 

(V-u h ,^)-(5'- 1 u h -w,^) + ((r + V-w)p ft ,( / 3 /l ) = (f,tp h ) for all ^6P (Tfe). (2.4) 

In what follows, we shall show an AMFEM algorithm based on the a posteriori error 
estimator developed in [21]. We note that our convergence analysis below is also valid for 
AMFEM based on the estimator proposed in [22] . 

Suppose that the module SOLVE outputs a pair of discrete solutions over Th, namely, 
(u/j,p/j) = SOLVE(7h). The estimator in [21] consists of several indicators with different 
weight factors, where the elementwise estimator f]^- h (uh,Ph, K) can, for convenience, be 
abbreviated to 

rf Th (u h ,p h ,K) :=Dlhl\\S-W\\ 2 K + o? K hl\\R K \\l+ £ ^llbt^" 1 ^)]^. 

aee K 

Here olk — mm(hK/y/cs,K, l/V'Vr^)- Rk is the elementwise residual defined by 
Rk ~ f - V • u h + (S^Uh) ■ w - (r + V • w)p h , 

and Dk, D a denote two variants of coefficients over each element K e Th and each side 
a G Eh respevtively given by 

D K := c w ^ K + C w r K a K + max c w . r ,K' + max , 

K':K'V\K^% K'-.K'nK^d) CwrK' 


1 1 r 2 h 2 r< 2 

u := - max Cs # + - mm{ max , max ). 

2 K:Xna-#0 2 KiifriCT^e Cvi, r ,K K-.Kn&^V) cs,k 


Define the global and local errors, and £k, of the stress and displacement variables as 


£%, £ 2 K--=\\S- 1/2 (n-u h )\\ 2 K + c w , r>K \\p- Ph \\ 2 K . (2.5) 

Ken 


From [21], or [22] but with different forms of Dk and D a , it holds an upper bound estimate 

e 2 h < Ctfl := C x i]\- h {\\ h ,p h ,Th) ■= Ci VT h ( u h>Ph,K), (2.6) 

Ken 

where the positive constant C\ depends only on the shape regularity of the meshes. 

For a given triangulation Th and a pair of corresponding discrete solutions (\ih,Ph) G 
RTo(Th) x Po(Th), we assume that the module ESTIMATE outputs the indicators 

{rf Th {u h ,p h ,K)} Kerh = ESTIMATE (u^ ,Ph,Th)- 
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Let Rk denote the mean of R K over each element K E Th- We define the oscillation 

osc 2 := h 2 K \\R K - Rk\\ 2 k . (2.7) 
Ken 

We note that throughout this paper the triangulation Th means a refinement of Th, and 
all notations with respect to the mesh Th are defined similarly. We shall also use the notation 
A < B to represent A < CB with C > a mesh-size independent, generic constant. 

In MARK step, by Dorfler marking we select the elements to mark according to the 
indicators, namely, given a grid Th with the set of indicators {i]j- H ( u h,Ph, K)}kgt h and 
marking parameter 9 e (0, 1], the module MARK outputs a subset of making elements, 
Mh C Th, with 

M H = MAM({^ fl (u H , pH,K)} KeTH ,T H ,0) 
satisfying Dorfler property 

VTh( u h,Ph,M h ) ■= ( ^r„( u H,PH,K)) 1/2 > 6t] Th (uh,Ph,Th)- 

KeM H 

In REFINE step, we suppose that the refinement rule, such as the longest edge bisection 
[33, 34] or the newest vertex bisection [37, 28, 29], is guaranteed to produce conforming and 
shape regular meshes. Given a fixed integer b > 1, a mesh Th, and a subset Mh C Th of 
marked elements, a conforming triangulation Th is output by 

T h = REFINE(Th,Mh), 

where all elements of A4h are at least bisected b times. Note that not only marked elements 
get refined but also additional elements are refined to recover the conformity of triangulations. 

We now describe the AMFEM algorithm. In doing so, we replace the subscript H by an 
iteration counter called k > 0. Let To be a uniform triangulation with a marking parameter 
9 £ (0, 1]. The basic loop of AMFEM is then given by the following iterations: 

AMFEM algorithm 
Set k = and iterate 

(1) (u fe ,p fe )-SOLVE(r fe ); 

(2) UKu^p^K)}^ = ESTIMATE(u fe ,p fe ,T fc ); 

(3) M k = MARK({r 1 l(u k ,p k ,K)} Ke T k ,Tk,9); 

(4) T fc+ i = REFINE (7fc, M k ); k = k + l. 

We note that the AMFEM algorithm is a standard one in which it employs only the error 
estimator {77^ (u k ,p k , K)}x£T k and needs neither marking the oscillation nor the interior 
node property. 

3. Estimate for L 2 projection of the displacement. This section is devoted to the 
estimation of | \QhP—Ph \\, where Qh is the L 2 -projection operator onto Po(Th)- The estimate 
is one key to the proof of convergence without the quasi-orthogonality available due to the 
convection term. It gives as well a posteriori error estimates for the L 2 — projection of the 
displacement variable (see remark 3.1). 

Consider the following auxiliary problem: 


V • (SV<f>) + V<£ • w - r<f> = Q h p - p h in ft, 

(J) — on dfl. 


(3.1) 
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It is well known that there exists a unique solution <ft £ Hq (O) to the problem (3.1) when the 
convection and reaction terms satisfy r + 1/2V • w > (Assumptions (7)1) and (DA)) withs 
the following regularity estimate: 

U\\m = U\\m { UT h )<\\QhP-Ph\\. (3.2) 
Moreover, if O is convex, S £ C 1,0 (O) implies the estimate 

U\\m {U T h )<\\QhP-Ph\\. (3.3) 

We emphasize that we only need an estimate on ||(/)|| ff 2( if ) for each K £ Th, i.e., the as- 
sumption on S could be weaken in the sense that only (3.3) is required. In [9] Carstensen 
gave an example which shows that when S is piecewise constant, cf) satisfies (3.3) but is not 
77 2 -regular. 

Set z := £ 77(div, fl), and denote <f>h the L 2 -projection of 4> onto Po{%), and II ^ 
the interpolation operator from 77 (div, fi) onto RT (Th) with the following estimate: 

Hfc-^z-IIfcz)!! < M*i (U 7i) forall ze77(div,n). (3.4) 

We refer to [2, 7, 25] for the detailed construction of such an interpolation operator H h and 
the approximation property. 

From (2.1) and (2.3), we obtain 

(Q hP - Ph ,V -Il h z) = {p-p h ,V-U h z) = (,S- 1 (u-u ft ),n /l z). (3.5) 
An integration by parts implies 

(S-^u - u fc ), z) = (S-^u - u h ), = -(V • (u - u h ), 0). (3.6) 

From (2.2) and (2.4) it follows 

(V-(u-u h ),0fc) = (S*- 1 (u-u h )-w,^)-((r + V-w)(p-p /l ),0 /l ) 

= (5- 1 (u-u /l )-w,0 /l )-(p-p /M (r + V-w)^) (3.7) 
= (S^u - u ^) ' w,^) - {QhP-Ph, (r + V • w)0 h ). 

Denote 7 := (Q h (V(j) ■ w),Q h p - p h ) - {rfih, QhP - Ph)- In view of the commuting 
property of the interpolation operator Uh, a combination of (3.5)-(3.7) yields 

\\QhP-Ph\\ 2 = (Q hP - Ph ,Q h V ■ z) + I = (Q hP - Ph ,V ■ Il h z) + I 
= (S-Hu - u h ),n h z -z) + (S-\u - u h ), z) + 7 

= (S-^u - u h ),n h z - z) - (V • (u - lift), <t> - <Ph) - (V • (u - u fc ), <j> h ) + I 

= (S-^u - u h ),n h z - z) - (V • (u - u h ),4> - <j> h ) - (S-^u - u h ) ■ w, <j) h ) 

+(V • w(f> h , QhP - Ph) + (Qh(V</> • w), Q h p - Ph) 

= (S-^u - u h ),U h z - z) - (V • (u - u h ),<j> - 4>h) - (S-^u - u h ) ■ w, 4> h - 4>) 

-(5 _1 (u- u h ) ■ w,4>) + (V • w4> h ,QhP- Ph) + (Qft(V^ • w),Q h p-ph)- 

(3.8) 

Recall the postprocessed technique developed by VohraMk in [39], where a postprocessed 
approximation p/, to the displacement p is constructed such that — SK^Ph\K = and 
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tjo J K Phdx = Ph\k for all K e %■ Then, from w e RT (Th), we have 

-{S-\u-u h )-w,4>)= J2 f V(p-p h )-w0 

Ken Jk 

Ken Jk 

V / -V<j>-w(p-ph)+ (p-ph)w-n<j>-(V-w4>,p-ph) 

T 7~i- Jk JdK 


Ken 1 


IK JdK 

= -(V</> • w,p-p h ) - (V • w<p,p-p h ) [Ph]w ■ n<j>. 

o J a 

Notice that it holds 

(Qh(V<t>-w),QhP-Ph) = {Qh{V(f>-™),p-ph) = {Qh{V<p-w),p-ph) (3.10) 

and 

(V -w(p h ,Q h p-p h ) = (V • w<p h ,p-p h ) = (V • w<p h ,p-p h ). (3.11) 

For convenience, denote 

h ■= (5 _1 (u- u h ),U h z- z), I 2 := -(V • (u - u h ), 4> - 4> h ), 

h ■= -(5 _1 (u-Uh) -w,(p h -(p), h ■= -(V(f>-w-Qh{V(/>- w),p-p h ), 

h ■= -(V -w(0- 4>h),p-Ph), h-^-Yl [Ph]wn<p. 
From (3.8)-(3.1 1) we arrive at 

4 


IIQ/^-^ll 2 = {V -w(t),p- ph) + (\7 -w(j) h7 Q h p- p h ) + I 6 

»=i 

4 

= - (V • w(0- <t>h),p-Ph) ~ (V • w<p h ,p-p h ) (3.12) 

i=l 

6 

+ (V • w^, Q/jp - p h ) + 7 6 = ^ Jj. 


i=l 


In what follows, we separately estimate Jj, i = 1, 6. 

Lemma 3.1. Denote \\h\\ L oo^ the maximum norm of the mesh-size function h with 
respect to Th, eh the error defined in (2.5). Then it holds 

h<\\h\\ L oo [Q) e h \\Q h p-p h \\. (3.13) 
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Proof. From (3.4) and (3.3) it follows 

h = ^2(S- 1 (u-u h ),Il h z-z) K 

Ken 

< ]T \\S-^ 2 (u-u h )\\ K \\Il h z-z\\ K 

Ken 

^ \\h\\L°°(n)ehMm(\jT h ) 

< ||/i|U~(n) e h \\Q h p-Ph\\- 


(3.14) 


□ 


Lemma 3.2. It holds 

h < {\\h\\ L oo^e h + osc h )\\Q h p-p h \\. (3.15) 

Proof. Notice that (2.2) can be equivalently written as: 

{V-u,<p) K -(S- 1 u-w,<p) K + ({r + V-w)p,<p) K = (f,if) K for all tp e L 2 (K), K e T h . 

(3.16) 

Meanwhile, the relation (2.4) can be equivalently written as: 

(V-u/,, (S^Ufc-w, ip)K+{{r+V-Yr)ph, <p)k = (/, <p)k for all ip e P (K),K & %■ 

(3.17) 

For arbitrary <p e L 2 (K), let <^ denote the mean of ip over K e 7fc. A combination of 
(3.16) and (3.17) yields 

(V • (u - u h ),(p) K = (V • u, </?) K - (V • Ufc, <p - <pk)k - (V • u h , <pk)k 
= (R K , (p - <Pk)k + (S'" 1 (u - u h ) ■ w, <^)k - ((r + V • w)(p - pk)^)k 
= (Rk - Rk,<P~ <Pk)k + (S _1 ( u - u h ) ■ w,<p>) K - (0 + V • ™)(p-Ph), i P)K 

< (||i?K-i?xlk + ||s- 1 (u-u/ l )lkl|w|Uo 0(K) + c W:r ,Klb-p /l |k)lklk. 

Here we recall that Rk is the mean of the elementwise residual Rk over each K £ Th- The 
above inequality indicates 

1 1 Y7 ( Ml (V ■ {u-U h ),<p) K 

V • (U- Uh)\\K = sup 


Then it follows 


V eL^(K).^0 \m\L^(K) (3.18) 
< \\R k -Rk\\k+£k. 


h = - Yj (V-(u-u h ),^-0 fc ) 

< ^ ||V-(u-u fc )|kMW|k ^ 19 > 

< nifellr«x, (n) e fc +oscfc)||V0||. 


The desired result (3.15) follows from (3.19) and (3.2). □ 
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Lemma 3.3. It holds 

J3<IWIi~ ( ii)e/,||Qfcp-pfc||. (3.20) 

Proof. By noticing 

^3 = - X! (S _1 ( u - u fc) -Vf,4>h -4>)k 

Ken 

< \\S- 1/2 tu-*h)\\Kh K \\V<f>\\K ( 3 - 21 > 
Ken 

< \\h\\ L °°(n)e h \\<t>\\m 

the desired result (3.20) follows from (3.21) and the regularity estimate (3.2). □ 
Lemma 3.4. It holds 

h < \\h\\ L ^ { n) e h \\QhP-Ph\\- (3-22) 

Proof. Recall a local efficiency estimate of Kk \\S~ x xyh\\K as following (see Lemma 7.3 
in [21, 22]): 

h K \\S- X M h \\K<£K. (3.23) 
By triangle inequality we obtain 

h = - ^2 (\7<j)-w -Q h (\7(f)-w),p-p h ) K 
Ken 

-Ph\\K + \\Ph -Ph\\K)K (3.24) 

Ken 

< h*c\\<t>\\H'{K)(\\p-Ph\\K + hK\\S- 1 Uh\\K). 

Ken 
From (3.23) it holds 

Wp-PhWK + hKWS-'uhWK <£ K . (3.25) 
According to (3.24) and (3.25), we arrive at 

h< Y h K\\<t>\\H*(K)£K, 

Ken 

which, together with the regularity estimate (3.3) of </>, implies the result (3.22). 
□ 

Lemma 3.5. It holds 

h < e-h\\h\\ L <~(Q.)\\QhP ~Ph\\- (3.26) 
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Proof. From (3.25) it follows 

h = - Y ( v ' w (0- 4>h),P~Ph)K 

Ken 

^ Y \\<t>-4>h\\K{\\p-Ph\\K + \\Ph-Ph\\K) 

Ken 

< hKWVtPWKiWp-PhWK + hKWS-WWK) 

Ken 

< 


eh\\ri\\ L °°(n)\m\Hi, 

which, together with (3.2), yields the estimate (3.26). □ 
Lemma 3.6. It holds 

l6<e h \\h\\^ {n) \\Q h p-p h \\. (3.27) 

Proof. For any a € e° h , let denote the mean of <fi onto a, i.e., <f> a := - — - [ <pds. 

I°1 Jo 

According to the continuity of the means of traces of the postprocessed scalar ph (see Lemma 
6.1 in [39]), and noticing w e RT n (Th), we have 

/ \p h ]w ■ n<j) = / \p h ]xv ■ n(<j> - <j) a ) 

Jo Jo (3.28) 

< WiPh}\\o\\<P'M\o. 

A sidewise Poincare inequality and trace theory imply 

\\</>-4>*\\o < hoWltA^Mo < h a \\<f>\\ HHUuJa) . (3.29) 

From trace theorem, generalized Friedrichs inequality (see (2.2) in [39]), and the post- 
processed technique, we have 

\\\Ph}\\o = HtPfc] - 7^7 / ]ph]ds\\ a < h^WVhPhlU < ^HS-^fcH^. (3.30) 
M Jo 

A combination of (3.28)-(3.30) yields 

f \p h ]w ■ n<f> < h a \\S-^\\ u X /2 \\<l>\\H'(U u „y (3.31) 

J o 

In light of the local shape regularity of element, the above estimate leads to 

h < y k\\s-w\\.x /2 m\h^) 


oee 

< 


o 


(3.32) 


z,°°(n) 

Ken 


The desired result (3.27) follows from (3.32), (3.23) and (3.3). □ 
We now give an estimate of 1 1 QuP — Ph 1 1 ■ 

THEOREM 3.7. Let eh, osc^, and ||/i||z,oo/q) denote the error of the stress and displace- 
ment variables given in (2.5), the oscillation of data given in (2.7) , and the maximum norm 
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of the mesh-size function, respectively, with respect to Th- Then there exits a positive constant 
C 2 only depending on the shape regularity ofTh, such that 

1/2 

\\QhP-Ph\\ <C 2 C D (\\h\\ ] / OO ( Sl ) e h + 0SC h)> ( 3 - 33 ) 

where Cd is one variant of coefficients. 

Proof The estimate (3.33) follows from a combination of (3.12), (3.13), (3.15), (3.20), 
(3.22), (3.26), and (3.27). □ 

REMARK 3.1. A combination of the two estimates (3.33) and (2.6) actually give a pos- 
teriori bound for Q h (p — p h ). Furthermore, following [7], we denote 

^k NC ■= Uh G L 2 (n) : q h \ K G P k (K),VK E %, V / Ph4>ds = O,V0 G TZ k (dK)} 

K JdK 

for k > 1, and let p* h ,Ph be respectively the interpolates in C\j NC of the interelement La- 
grangian multiplier A/j and the displacement p ([7], pages 186-187; We note that in [7] u 
represents the displacement variable and p the stress variable). Following the same line as in 
[7], it holds the estimate 


\Ph 


-Ph\\ <\\hS- 1/2 (u-u h )\\ + \\Q h (p- Ph ) 


which gives an a posteriori error estimate for ph — p* h . 

REMARK 3.2. For a pure diffusion problem, i.e., w = r = in (1.1), it holds Ii — 0, i — 
3, 4, 5, 6. From the estimates of Ii and I 2 , we can obtain 

\\QhP-Ph\\<\\hS- 1 ^(u-u h )\\+ osc h , 

which results in the quasi-orthogonality 

(S-\n-u h ),u h - vl h ) < (ll^ 1/2 (u - u h )|| + osc h )||/ fc - f H \\, 

where we have used the fact that V • = Qhf = fh an d V • uh = Qnf = Jh- This 
estimate is somewhat different from the quasi-orthogonality results in [10, 13, 14, 19, 20]. 

4. Estimator reduction. Let uj a denote the patch of a e £/,, and define c w<T , D 2 K , 
Dj- h (K), Dj- h respectively by 

{maxfco 1 / 2 ,^, 1 / 2 ) if a = K D L, 2 ~ 2 

Cg 1 ^ if a e e K n dn, ^e K 

D\{K) := max(44c s i,C OK 4,^), D\ := £ D\(K), 

where Cdk in Dj- (K) is given by 


Ken 


C DK ■= 2max((C 4 C s / K + w \\ l °°(k)) , )• 

LEMMA 4. 1 . (Estimator reduction) For a triangulation Th with Mh C Th, let Th be 
a refinement of Th obtained by the AMFEM algorithm. Denote by D 2 ^ one variant of the 
coefficients onto the initial mesh To, and denote A := 1 — 2~ b / 2 , 

E 2 H -= E \\S- 1/2 (n h -UH)\\ K + c w ^ K \\ Ph -pH\\ 2 K . (4.1) 

Ken 
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Then for any < 5 < 1, it holds 

(4.2) 

Proof. By triangle inequality and Young inequality, we have 
^llS-^ll^^a + JJ^IIS-^ffll^ + a + J-^^c^llS-Va^-u^l^. (4.3) 
Inverse inequality implies 

||V • (u h - u h )\\k < CiCy^h^WS- 1 ' 2 ^ - u H )\\ K , 

which leads to 

h 2 K \\f - V • u h + S'W ■ w - (r + V • w)p h \\ 2 K 

< (1 + S)h 2 K \\f -V-u H + S^uh ■ w - (r + V • w)p H \\ 2 K (4.4) 
+ (1 + tf-^CDArdlS-Va^ - u H )\\ 2 K + c w ^ K \\ Ph - p H \\ 2 K ). 

From inverse inequality and local shape regularity of the mesh, it follows 
which, together with Young inequality, yields 

Mfrt^Ufc)]^ <(l+^)^||[7t.(5- 1 u ff )]||2 + (l + 5-i)C|4j|5- 1 / 2 (u, l -u ff )||^. 

(4.6) 

Summing (4.3) and (4.4) over all elements K e 7/j, summing (4.6) over all interior sides 
a G e°, and noticing the monotonicity of variants of coefficients, we obtain 

rj?r h (u h ,p h ,T h ) < (1 + ^^(ujr.pH.Tfc) + C 2 (l + S^D^E 2 ,. (4.7) 

For a marked element K e M H , we set Tu.k ■= {K 1 e Th|if' C X}. It holds 

f ^ V 2 - h (u H ,p H ,K')<2- b / 2 r ] 2 - H (vi H ,p H ,K) for K eT H /T h , 
{ K'en, K 

[ rfi- h (u H ,PH,K) < Vt h ( u h,Ph,K) for KeT H /M H , 

which results in the following estimate 

V?r h {uH,PH,T h ) < 2- b / 2 V ^ H (u H ,p H ,M H )+VT H ( u H,PH,TH/T h ) 

(4.8) 

= Vt h ( u H'Ph,Th) - \ii^ H (u H ,p H ,M H )- 
The desired result (4.2) follows from (4.7), (4.8) and the monotonicity of Dj- h . 

5. Proof of theorem 1.1. . In this section, we show that the error plus some quantity 
uniformly reduces with a fixed factor on two successive meshes, which shows AMFEM is 
convergent. 

LEMMA 5.1. Let {u-hiPh) € RTn{Th) x PoiTh) be the approximation solutions to 
the stress and displacement variables with respect to 7~h, and eh the error of the stress and 
displacement variables with respect to Th- Denote by h the mesh-size of the quasi-uniform 
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initial mesh %, by D\ one variant of the coefficients determined by Co- Then it holds for 
any < 5\ < 1 

HS-V^u - U/OH 2 < 25^D 1 (h a e 2 h + h 2 (\\f - f h \\ 2 + HW^S" 1 ^ • w)|| 2 )) 

+ ||S'- 1 /2( U _ Uff) ||2 _ ||5-V2( Ufc - u H )\\ 2 + 1/25!||V • (u fc - u H )|| 2 . 

(5.1) 

Proof. From (2.1) and (2.3) we get 

HS-V^u-OH 2 = ||^- 1 /2(u-u if )|| 2 -||5- 1 / 2 K-u ff )|| 2 
-2(,S- 1 (u-u, l ),u h -u fl ) 
= \\S-V*(u-u H )\\*-\\S-V*(u h -u H )\\* 
-2(p-p/,,V- (u/, -u ff )) 

= ^-^(u-u^lp-H^K-u^)!! 2 

-2(Q h p -p h ,V ■ (u h - u H )). 
The definition of the residual Rk and the assumptions of data imply that for each K e Th, 

h K \\R K - Rk\\k < CeihlWViS-W ■ w)\\ K + h K \\f - f h \\ K ), 

which, together with the fact ||/i||l~(^) < h and the definition (2.7) of osc h , yields the 
estimate 

osc 2 < 2Clh 2 (\\f - M| 2 + WhVhiS-W ■ w)|| 2 ). (5.3) 

Applying the above estimate (5.3) to (3.33), we obtain 

\\QhP-Ph\\ 2 < Dtihoek + h 2 (\\f - M| 2 + WhWniS-W ■ w)|| 2 )). (5.4) 

In light of Young inequality, we have 

2(Q hP - Ph ,V- (u fc - u H )) < 2S^\\Q hP - Ph \\ 2 + ^/2||V • (u h - u H )\\ 2 . (5.5) 

The desired result (5.1) follows from a combination of (5.2), (5.4), and (5.5). □ 

LEMMA 5.2. Let D 2 and D3 be two variants of the coefficients respectively given by 

D 2 := max \\vr\\ 2 L oo {K) Cg* K , D 3 := max G^^ K . 

Under the assumptions of Lemma 5.1, it holds 

||V • (u- u,0|| 2 < ADMhvel + h 2 (\\f - M| 2 + WhVhiS-W ■ w)|| 2 )) 
+ ||V • (u - u ff )|| 2 - 1/2||V • (u fc - u H )|| 2 + 4C 2 ||5- 1 / 2 (u- U/l )|| 2 . 

(5.6) 

Proof. Notice 

||V-(u-0|| 2 - ||V-(u-u if )|| 2 -||V-(u ft -u H )|| 2 -2(V-(u-u h ),V-(u, l -u ff )). (5.7) 
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A combination of (2.2) and (2.4) leads to 

(V • (u - Ufc ), V • (u fc - u H )) = (S-^u - Ufc ) ■ w, V ■ (u h - u ff )) 

-((r + V ■ w)(g hP - p h ), V • (ufc - u H )) (5.8) 

< 2C 2 ||5- 1 /2( U - U/l )|| 2 + 2D 3 ||Q fc p-p h || 2 + 1/4||V • (u fc - u H )\\. 

The estimate (5.6) follows from (5.7), (5.4) and (5.8). □ 

LEMMA 5.3. Let be one variant of the coefficients given by D4 :— m&XKeT h c w.r,K- 
Under the assumption of Lemma 5.1, it holds 

^2 Cw,r,K\\V-Ph\\ 2 < ^2 Cw^kWp-PhWk - 1/2 X] C ™,r,K\\Ph- PH\\ 2 K 

KeT h KeT H KeT H 

+ 2D A D x {h Q e\ + h 2 (\\f - M| 2 + HWfc^- 1 ^ • w)|| 2 )). 

(5.9) 

Proof. For each element K e Th, it holds the following identity 

\\p-Ph\\ 2 K = \\p-Ph\\k -Wph-PhWk -2(p-Ph,Ph-PH) K 

= \\p-Ph\\ 2 k - \\Ph -Ph\\ 2 k - 2(QhP-Ph,Ph-PH)K- 

Notice that c W: , jf does not change from Th to Th- Summing (5.10) by multiplying c w>r ^ 
over all elements ]f e Tk, we have 

C ™,r,K\\P ~ Ph\?K < ^^KWP - Ph\\k - 1/2 51 C w,r,/f|K -P//Hk 

ifeTh atgTh xeTff 

+ 2D 4 \\Q hP - Ph \\ 2 . 

(5.11) 

The conclusion (5.9) follows from (5.11) and (5.4). □ 

In what follows, we show the reduction of the error. To this end, set 71,72, £0, an d ^1 to 
be any positive constants, which will be determined below. Introduce the following quantity: 

A 2 := ( 5 1 (l- £o )- 1 ||V.(u-u, l )|| 2 +7i^+72(||/-M| 2 + ||W, l ( > S- 1 u,.w)|| 2 ), (5.12) 

where fh is the i 2 — projection of / onto Po(Th). We note that the definition of Ah is similar 
to A h . 

Theorem 5.4. Let(u h ,p h ) e RT (T h ) x P (T h ) and{u H ,p H ) G RT {T H ) x P {T H ) 
be the approximation solutions to the stress and displacement variables with respect to Th and 
Th, respectively. Denote by eh and en the errors of the stress and displacement variables 
with respect to Th and Th, respectively. Let ho be the mesh-size of the quasi-uniform initial 
mesh To, and q and a G (0, 1) two constants to be determined below. Then it holds 

e 2 + (1 - h q)A 2 h < a 2 (e 2 H + (1 - h q)A 2 H ) (5.13) 

when h < \^ 1 

Proof. For convenience, denote 

D 5 := 2D 1 + A5lD x Dz + 2<5 1 £> 1 £> 4 . 

Recalling the definition, (4.1), of E H , a combination of (5.1), (5.6) and (5.9) indicates 

e 2 + <5i||V • (u - u h )\\ 2 < e 2 H + 5i||V • (u - u ff )|| 2 - 1/2E% + h D 5 S^e 2 

+A5 1 D 2 e 2 h + D 5 6^h 2 (\\f - f h \\ 2 + HW^V • w)|| 2 ). 

(5.14) 
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For a constant e > which will be determined below, denote e := j^. We firstly choose 
Si with A5iD 2 < s . The reliable estimate, (2.6), of the stress and displacement variables, 
together with (5.14), implies 

e2 + 5i||V-(u-u fc )|| 2 < e 2 H + 5 1 \\V-(u-u H )\\ 2 -l/2E 2 H + C 1 D 5 5^ 1 h r ] 2 h 

+e el + D 5 S^h 2 (\\f - M| 2 + HWfc^- 1 ^ • w)|| 2 ), 

which results in the following inequality: 

4 + T^IIV • (u - u,)|| 2 < (1 + e)e\ + • (u - u ff )|| 2 - 

+m^r) h ovl + Mi^^dl/ - Ml 2 + \\hv h (s-w ■ w)|| 2 ). 

_ (5.15) 

According to triangle inequality and inverse inequality, it holds, for each K e Th, 
hKWViS^Uh ■ w)\\ K < hKWViS^UH ■ w)\\ K + CVIIS" 1 ^ - u H ) ■ w\\ K . (5.16) 

Notice that || / - f h \\ K < \\f - f H \\K for all K e Th- For any given <5 3 > which will be 
determined below, (5.16) and Young inequality imply 

11/ - Ml 2 + WhVhiS-W ■ w)|| 2 < (1 + S 3 )(\\f - IhW 2 + WHVniS-'uH ■ w)|| 2 ) 

+ (1 + 5^ 1 )C^D 2 \\S- 1 / 2 (u h - u H )\\ 2 . 

(5.17) 

From the definition, (5.12), of A 2 h , the estimator reduction (4.2) with the marking strat- 
egy, the estimates (5.15) and (5.17), and the fact \\S^ 1/2 {u h - u H )\\ 2 < E%, it holds, for 
any given 5 2 > which will be determined below, 

el + A\ < (1 + e)e 2 H + • (u - u ff )|| 2 - ^ife + s^yM 2 . 

+C7 3 2 (1 + S 2 1 )D 2 Toll E 2 H + 7l (l + S 2 )(1 - X0 2 ) V 2 H+l2 C 2 D 2 (l + j- 3 )E 2 H 
+ MT^(H/ " Ml 2 + WhVniS-W ■ w)|| 2 ) 
+72(1 + 5 3 )(\\f - / H || 2 + WHVHiS-'uH • w)|| 2 ). 

(5.18) 

We next choose 71 and 72 such that 

7^(1 + S 2 ^ To = , 2 C 2 D 2 (1 + ^ = 

Then it follows 

el + A 2 h < (l + e)e 2 f + T ^-||V.(u-u ff )|| 2 +7i(l + <5 2 )(l-A0 2 )r ? 2 f 

+ MT^o(ll/ - Ml 2 + WhVniS-W • w)|| 2 ) + s$^hovt 
+72(1 + S 3 )(\\f - .fell 2 + WHVniS^UH ■ w)|| 2 ). 

(5.19) 

For any given 84,65 > which will be determined below, the reliable estimate (2.6) on 
Th, i-e. e 2 H < C\r^ H , and the above estimate (5.19), indicate 

e 2 h + A 2 < (1 + e - l/2A0 2 7l (l + 5 2 )C^)e 2 H + 7l (l + S 2 )(l - l/2Xe 2 ) V 2 H 

+ M^M|| V . (u _ Uff) || 2 + j-^^dl/ - Ml 2 + \\hV h (S-W ■ w)|| 2 ) 

+ I^%h oV 2 + (1 - <5 5 ) 72 (1 + £ 3 )(||/ - f H \\ 2 + \\HV H (S-'u H ■ w)|| 2 ) 

+ ^||V • (u- u ff )|| 2 + (1 + ^572(11/ - .fell 2 + WHVniS-'uH ■ w)|| 2 ). 

(5.20) 
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Now we fix a sufficiently small 82 and, subsequently, a sufficiently small e such that 

C \0 2 
a 2 := max(l + e - -^Afl 2 7 i(l + 8 2 ), (1 + «5 2 )(1 - — ), 1 - 8 4 , (1 - <J 6 )(1 + S 3 )) < 1. 

Let D 6 be one variant of the coefficients given by 

A :=4 naeKdlwIlioc^Cg^^^^c-^,!). 

From (3.18), we get 

||V ■ (u - u ff )|| 2 < D 6 (\\f - f H \\ 2 + WHVniS-'uH • w)|| 2 + r, 2 H ). (5.21) 
We further choose Si(i = 3, 4, 5) such that 

t— — — < 72^0, (1 + W < C 8 h . 

1 — £0 

In fact, we may firstly fix 5% satisfying <5 3 < C$h , then choose 65 such that 

TT^ mm(1 ' T+aa >" 

Finally, by noticing the choice of 72, we can choose <5 4 with 

j. . , fto 1 - £0 1 \ 
04 < mm 1, =— ). 

A <*i 4C 2 ^ 2 (l- £o )(l + ^ 1 ) 
These choices, together with (5.20) and (5.21), lead to 

el + A 2 < a 2 (e 2 H + A 2 H ) + J ^h oV 2 

+ sj^o) h o(\\f Ml 2 + \\hV h (S-'u h ■ w)|| 2 ) 
+72/10(1 + Cs)(\\.f - .fe|| 2 + WHWHiS^UH ■ w)|| 2 ) + j 2 h oV %. 

Let q be one variant of the coefficients given by 

(n , n 72 giA 1 A ftp. 

g := max C 8 + 1, — , — r — , -z— r— . 

71 di(l-e )7i oi(l-£o)72 

From (5.22) we arrive at 

e 2 h + A 2 < a 2 (e 2 H + A 2 H ) + qh (A 2 + A 2 H ), 

which implies 


(5.22) 


e 


2 h + (l- qh )A 2 < a 2 e 2 H + (a 2 + qh )A 2 H . (5.23) 


We finally choose ho such that 

a 2 + qh 2 1 + & 2 

0<-j 7— < a := — o — , 

1 — qho 2 

which yields the assertion (5.13) with h < (1 - a 2 )/(g(l + a 2 )). □ 
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REMARK 5.1. (Choices of the initial mesh size) Some simple calculations show 

/ r n/j x \ D 5 h C 7 2 

q < m&x{D(5 1 ,6 2 ), -jr- } 

4di(l - s yD 2 


with 


C 2 (l + S^D 2 - C 2 D<, 
D(8 U 5 2 ) := max{C 8 + 1, 3( 2 ' T \ — ^ 5 }. 


Then it holds 


which indicates 


1 . , 1 AS 1 (l-e ) 2 C 2 D 2 

^ mm{ )3(W } ' 


1 1 , , 4S 1 (l-e ) 2 C^D 2 

q D{5 U 8 2 ) D b 

required in Theorem 5.4, the initial mesh size h is assumed to satisfy h < 
Then eventually we may choose h with 

, , . s I -a 2 45 1 (l-e ) 2 C 2 D 2 . 
h < mm{— — }. 

The proof of Theorem 1.1 Theorem 5.4 shows that the error of the stress and displace- 
ment variables plus the quantity A 2 , uniformly reduces with a fixed factor a 2 between two 
successive meshes. Replace the subscripts H and h respectively by the iteration counters k 
and k + 1, we then obtain Theorem 1.1 directly from Theorem 5.4. 

6. Numerical experiments. In this section, we test the performance of the adaptive 
algorithm AMFEM described in section 2 with four model problems. We are thus able to 
study how meshes adapt to various effects from lack of regularity of solutions and convexity 
of domains to data smoothness, boundary layers and changing boundary conditions. We note 
that the implementation of AMFEM is done without enforcing the interior node property in 
the refinement step. 

6.1. Model problem with singularity at the origin. We consider the problem (1.1) in 
an L-shape domain fl = {(-1,1) x (0,1)} U {(-1,0) x (-1,0)} with w = r = and 
/ = 0. The exact solution is given by 

p(p,0) =p 2/3 sin(2#/3), 

where p, 8 are the polar coordinates. 

Since this model possesses singularity at the origin, we see in the left figure of Fig 6.1 
that the refinement concentrates around the origin, which means the predicted error estimator 
captures well the singularity of the solution. The right graph of Fig 6. 1 reports the estimated 
and actual errors of the numerical solutions on uniformly and adaptively refined meshes. It 
can be seen that the error of the stress and displacement in L 2 norm uniformly reduces with 
a fixed factor on two successive meshes after several steps of iterations, and that the error on 
the adaptively refined meshes decreases more rapidly than the one on the uniformly refined 
meshes. This shows that the adaptive mixed finite element method is convergent with respect 
to the energy error. 
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FIG 6. 1 . A mesh with 14692 triangles (left) and the estimated and actual errors in uniformly / adop- 
tively refined meshes ( right) for the marking parameter 9 = 0.5. 


6.2. Model problem with inhomogeneous diffusion tensor. We consider the problem 
(1.1) in a square domain = ( — 1, 1) x (—1, 1) with w = r = and / = 0, where fi is 
divided into four subdomains 17,; (i = 1, 2, 3, 4) corresponding to the axis quadrants (in the 
counterclockwise direction), and the diffusion-dispersion tensor S is piecewise constant with 
S = Sil in fii. This model problem is taken from [24, 36, 39]. We suppose the exact solution 
of this model has the form 

p(p,9) — p a (disin(a0) + biCos(aO)) 

in each fij with Dirichlet boundary conditions. Here p, 9 are the polar coordinates in SI, 
cij and bi are constants depending on fij, and a is a parameter. We note that The stress 
solution, u = — iSVp, is not continuous across the interfaces, and only its normal component 
is continuous. It finally exhibits a strong singularity at the origin. We consider two sets of 
coefficients in the following table: 


Case 1 

Case 2 

Si — S3 — 5, S 2 — S 4 — 1 

Si — S3 — 100, S2 — S4 — 1 

a = 0.53544095 

a = 0.12690207 

ai = 0.44721360, bi = 1.00000000 
a 2 = -0.74535599, b 2 = 2.33333333 
a 3 = -0.94411759, b 3 = 0.55555555 
a 4 = -2.40170264, b 4 = -0.48148148 

ai = 0.10000000, fci = 1.00000000 
a 2 = -9.60396040, b 2 = 2.96039604 
a 3 = -0.48035487, b 3 = -0.88275659 
a 4 = 7.70156488, b 4 = -6.45646175 


In MARK step, the marking parameter 9, in terms of Dorfler marking, is chosen as 0.7 
in the first case and as 0.94 in the second case. Table 6.1 shows for Case 1 some results of 
the actual error e%, the a posteriori indicator rjk, the experimental convergence rate, EOC^, 
of Ek, and the experimental convergence rate, EOC^, of i]k, where 

Fnr log(e fc _i/e fc ) logfa-i/V) 

log(DOF fc /DOF fc _i)' £jU °"' logpOFfc/DOF^)' 

and DOF/j denotes the number of elements with respect to the k— th iteration. We can see 
that the convergence rates EOC# and EOC ?) are close to 0.5 as the iteration number k = 15, 
which means the optimal decays of the actual error and the a posteriori error indicator r)k 
are almost attained after 15 iterations with optimal meshes. 

Fig 6.2 shows an adaptively refined mesh with 4763 elements and the estimated and 
actual errors against the number of elements in adaptively refined meshes for Case 1. Fig 6.3 
shows an adaptively refined mesh with 1093 elements and the actual error against the number 
of elements in adaptively refined meshes for Case 2. 
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Table 6.1 

Results of actual error E^, a posteriori indicator rjk, and their convergence rates EOC e and EOCrj - Case 1 


k 

DOF fc 

e-h 

Vk 

EOCb 

EOC„ 

1 

8 

1.3665 

5.0938 



2 

20 

1.1346 

3.4700 

0.2030 

0.4189 

9 

2235 

0.1776 

1.1115 

0.4016 

0.4004 

11 

7165 

0.1106 

0.7111 

0.3851 

0.4004 

12 

13188 

0.0871 

0.5566 

0.3915 

0.4015 

14 

43785 

0.0510 

0.3365 

0.4707 

0.4476 

15 

76770 

0.0387 

0.2581 

0.4915 

0.4724 



From the left figures of Fig 6.2-6.3, we can see that the refinement concentrates around 
the origin, which means the AMFEM algorithm detects the region of rapid variation. In 
the right graphs of Fig 6.2-6.3 each includes an optimal convergence line, which shows in 
both cases, the energy error performs a trend of descending with an optimal order convergent 
rate after several steps of adaptive iterations for the problem with strongly discontinuous 
coefficients. We note that the energy error is approximated with a 7-point quadrature formula 
in each triangle. 

6.3. Convection-dominated model problem: boundary layer. In this example, we 
take f2 = (0, 1) x (0, 1) in W 2 , and choose w = (1, 1) and r = 0. Further, we set p = on 
dft, and select the right-hand side / such that the analytical solution to (1.1) is given by 

exp(-f)-l exp(-f)-l 

The solution is smooth, but has boundary layers at x = 1 and y = 1, with layer width of 
order 0(e). This problem is well-suited to test whether the estimator is able to pick up the 
steep gradients near these boundaries. 

We start computations from the origin mesh consisted of 8 right-angled triangles, and we 
choose the marking parameter 9 — 0.5 in the adaptive algorithm AMFEM. 

Fig 6.4 shows the mesh with 10838 triangles (left) and the postprocessing approximation 
to the scalar displacement p on the corresponding adaptively refined mesh (right) in the case 
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FIG 6.3. A mesh with 1093 triangles (left) and the actual error against the number of elements in 
adoptively refined mesh ( right): Case 2. 



FIG 6.4. A mesh with 10838 triangles (left) and postprocessing approximate displacement on the 
corresponding adoptively refined mesh (right) for e — 0.01. 


e = 0.01. Here the value of the postprocessing approximation on each vertex is taken as the 
algorithmic mean of the values of the displacement finite element solution on all the elements 
sharing the vertex. The reason for the postprocessing is that the displacement finite element 
solution is not continuous on each vertex of the triangulation. We see that the refinement 
focuses around boundary layers, which indicates that the estimators actually capture boundary 
layers and resolve them in convection-domianed regions. In addition, the postprocessing 
approximation to the scalar displacement obtains satisfactory results. 

Fig 6.5 shows the actual error (energy error) results against the number of elements in 
adaptively refined meshes for e = 0.1 (left) and e — 0.01 (right), including two theoretically- 
optimal order (-1/2) convergence lines. We see that in each case the actual error descends 
almost at the optimal rate of convergence after several steps of iterations. The numerical 
results confirm our theoretical analysis. 

6.4. Convection-dominated model problem: interior and boundary layer. Set the 

domain O = [—1,1] x [—1,1] with non-homogeneous Dirichlet boundary conditions, the 
velocity field w = (2, 1), and the reaction term r = in (1.1). The source term / = 0, the 
Dirichlet boundary conditions are as follows: p = along the left and top sides of the square 
andp =100 along the right and bottom sides. The exact solution of this problem is unknown, 
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Number of triangles Number of triangles 


FIG 6.5. Actual error against the number of elements in adaptively refined meshes for e = 0.1 (left) 
and e = 0.01 f right) for the marking parameter 8 — 0.5. 


but it is known that it exhibits an exponential boundary layer at the boundary x = 1, y > 
and a parabolic interior layer along the line connecting the points (—1,-1) and (1,0). 



FIG 6.6. A mesh with 47324 elements (left) for e = 0.01, 6 = 0.8 and estimated error against the 
number of elements in adaptively refined meshes (right) for s = 0.1, 9 = 0.5. 

We still perform the AMFEM algorithm described in section 2 from the origin mesh con- 
sisted of 8 right-angled triangles. From the left graph of Fig 6.6, we can see that when using 
adaptive refinement the mesh concentrates close to the exponential and parabolic layers. We 
note that the refinement first occurs close to the region x = 1, y > 0, since the exponential 
layer is more stronger than the parabolic layer. The left graph also illustrates that the a pos- 
teriori error estimator exactly capture the behavior of the solution. The right graph of Fig 6.6 
shows that the estimated error rapidly reduces starting from the fourth step of iterations, and 
reaches the optimal rate (-1/2) of convergence until the seventeenth step. This convergence 
result is consistent with our theoretical analysis. 

REFERENCES 

[1] M. Anisworth and J.T.Oden, A posteriori error estimation in finite element analysis, Wiley, Chichester, 2000. 

[2] D.N. Arnold, R.S.Falk and R.Winther. Finite element exterior calculus, homological techniques, and applica- 
tions. Acta Numerica, pages 1-155, 2006. 

[3] I. Babuska and W. Rheinboldt, A posteriori error estimates for the finite element method, Internat. J. Numer. 
Methods Engrg, 1978, 12, 1597-1615. 


22 


[4] I. Babuska and W. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. 
Anal, 1978, 15, 736-754. 

[5] I. Babuska and T. Strouboulis, The finite element method and its reliability, Clarendon Press, Oxford, 2001. 
[6] R. Becker, S. Mao, An optimally convergent adaptive mixed finite element method, Numer. Math, 2008, 111: 
35-54. 

[7] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer, Berlin-Heidelberg-New York, 
1991. 

[8] J.M. Cascon, C. Kreuzer, R.H. Nochetto and K.G. Siebert, Quasi-optimal convergence rate for an adaptive 

finite element method, SIAM J. Numer. Anal, 2008, 46(5): 2524-2550. 
[9] C. Carstensen, A posteriori error estimate for the mixed finite method, Math. Comp., 66(218), 1997, 465-476. 
[10] C. Carstensen and R.H.W. Hoppe, Error reduction and convergence for an adaptive mixed finite element 

method, Math. Comp, 2006, 75, 1033-1042. 
[11] C. Carstensen and J. Hu, A unifying theory of a posteriori error control for nonconforming finite element 

methods, Numer. Math., 2007, 107, 473-502. 
[12] C. Carstensen, J. Hu, A. Orlando, Framework for the a posteriori error analysis of nonconforming finite 

elements, SIAM J. Numer. Anal., 2007, 45, 6882. 
[13] C. Carstensen and H. Rabus, An optimal adaptive mixed finite element method, Math. Comp., 80 (2011), 

649667. 

[14] L. Chen, M. Hoist and J.C. Xu, Convergence and optimality of adaptive mixed finite element methods, Math. 

Comp., 2009, 78(265) :35-53. 
[15] R G. Ciarlet, The finite element method for elliptic problems. Nort-Holland, Amsterdam, 1978. 
[16] A. Demlow and R. Stevenson, Convergence and quasi-optimality of an adaptive finite element method for 

controlling I? errors, Math. Comp., 2011, 117 :185-218. 
[17] J. R. Douglas and J. E. Roberts, Global estimates for mixed methods for secod elliptic equations, Math. 

Comp., 44 (1985), 39-52. 

[18] S.H. Du and X.R Xie, Residual-based a posteriori error estimates of nonconforming finite element method 
for elliptic problem with Dirac delta source terms, Science in China Series A: Mathematics, 2008, 51(8): 
1440-1460. 

[19] S.H. Du and X.R Xie, Error reduction, convergence and optimality of an adaptive mixed finite element 
method, Journal of systems science and complexity, 201 1, 24: 1-14. 

[20] S.H. Du and X.R Xie, Error reduction, convergence and optimality for adaptive mixed finite element methods 
for diffusion equations, Journal of Computation Mathematics, 2012, 30 (5): 483-503. 

[21] S.H. Du and X.R Xie,A new residual-based posteriori error estimators for lowest-order Raviart-Thomas ele- 
ment approxiamtion to convection-diffusion-reaction equations, J. Comput. Math., revision resubmitted 

[22] S.H. Du, A new residual posteriori error estimates of mixed finite element methods for convection-diffusion- 
reaction equations, Numerical methods of PDE, accepted. 

[23] W. Dorfler, A convergent adaptive algorithm for Poisson's equation, SIAM J. Numer. Anal,1996, 33, 1106- 
1124. 

[24] G. T. Eigestad and R. A. Klausen, On the convergence of the multi-point flux 0-mefhod: Numerical ex- 
periments for discontinuous permeability, Numer. Methods Partial Differential Equations, 21 (2005), 
1079-1098. 

[25] R. Hiptmair, Canonical construction of finite elements. Mathematics of Computation, 1999, 68: 1325-1346. 
[26] J. Hu, Z.C. Shi and J.C. Xu, Convergence and optimality of the adaptive finite element method. Research 

Report, 19(2009), School of Mathematical Science and Institute of Mathematics, Peking University, 

available at www.math.pku. edu.cn:800/var/preprint/7197.pdf; Numer. Math., to appear. 
[27] K. Mekchay and R.H. Nochetto, Convergence of adaptive finite elemet methods for general second linear 

elliptic PDEs, SIAM J.Numer. Anal, 2005, 43: 1803-1827. 
[28] W.F. Mitchell, A comparison of adaptive refinement techniques for elliptic problems. ACM Trans. Math. Soft, 

1989, 15(4): 326-347. 

[29] W.F. Mitchell, Optimal multilevel iterative methods for adaptive grids. SIAM Journal on Scientific and Sta- 
tistical Computing, 1992, 13: 146-167. 

[30] P.Morin, R.H. Nochetto and K.G. Siebert, Data oscillation and convergence of adaptive FEM, SIAM J. Numer. 
Anal, 2000, 38(2): 466-488. 

[31] P.Morin, R.H. Nochetto and K.G. Siebert, Local problems on stars: a posteriori error estimators, convergence, 

and performance, Math. Comp, 2003, 72 (243): 1067-1097. 
[32] P. Morin, R.H Nochetto and K.G. Siebert, Convergence of adaptive finite element methods, SIAM Review, 

2002, 44(4), 631-658. 

[33] M.C. Rivara, Mesh refinement processes based on the generalized bisection of simplices. SIAM J. Numer. 
Anal, 1984, 21: 604-613. 

[34] M.C. Rivara, Design and data structure for fully adaptive, multigrid finite element software. ACM Trans. 

Math. Soft, 1984, 10 :242-264 
[35] P. A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathe- 


CONVERGENCE OF ADAPTIVE MIXED FINITE ELEMENT METHOD 


23 


matical Aspects of Finite Element Methods (Proceedings of the Conference of the C.N.R., Rome, 1975), 
Lecture Notes in Math. 606, Springer, Berlin, 1977, 292-315. 

[36] B. Riviere and M. F. Wheeler, A posteriori error estimates for a discontinuous Galerkin method applied to 
elliptic problems, Comput. Math. Appl., 46 (2003), 141C163. 

[37] E.G. Sewell, Automatic generation of triangulations for piecewise polynomial approximation. In Ph. D. dis- 
sertation. Purdue Univ., West Lafayette, Ind., 1972. 

[38] R. Verfiirth, A review of a posteriori error estimation and adaptive mesh-refinement techniques, Wiley- 
Teubner, 1996. 

[39] M. Vohrah'k, A posteriori error estimates for lowest-order mixed finite element discretizations of convection- 
diffusion-reaction equations, SIAM J. Numer. Anal., 45(4), 2007, 1570-1599. 


